Analytical solutions for cavity contraction in strain-softening materials with linear or exponential strength decay

This paper presents an analytical investigation into the contraction of spherical and cylindrical cavities excavated in strain-softening rock masses obeying the Mohr–Coulomb or Tresca yield criterion, with linear or exponential uniaxial compressive strength decay. The derivation of the ground response curves is based on the simplifying assumption that the strains inside the plastic zone are completely plastic. This significantly simplifies the mathematical formulation, enabling the derivation of closed-form solutions. An alternative simplifying approach which partially neglects the elastic strains inside the plastic zone and which is commonly adopted in the literature, is also examined. The accuracy of the simplified solutions is evaluated by comparing their predictions with rigorous solutions obtained by numerical finite-difference analyses. The investigation demonstrates that the proposed closed-form solutions represent a significant improvement on those based on the commonly-made simplifying assumption involving partial neglect of elastic strains.


List of symbols
The contraction of cavities excavated in strain-softening rock masses has so far been studied using analytical, semi-analytical and numerical methods, which can generally be categorized into two groups based on their treatment of strain within the plastic zone: rigorous methods utilizing incremental plasticity and simplified methods.
Recent works and comprehensive reviews of the state of the art can be found in [1][2][3][4][5][6][7][8][9][10][11] .From a theoretical perspective, the rigorous methods are robust, but expressing closed-form solutions for the cavity contraction problem has proven challenging.Instead, numerical integration is typically required.In practice, for arbitrary yield conditions and flow rules, the system of the governing ordinary differential equations can only be solved numerically, using for example the Runge-Kutta algorithm [12][13][14] .In general, these methods are too complex for practical use, due to the need to solve intricate second-order (or systems of) differential equations or to implement advanced, tailor-made constitutive models in existing finite element codes.Closed-form solutions are by far easier to apply in engineering design computations, but exist only for specific constitutive models and boundary value problems as can be seen from Table 1, which provides an overview of existing analytical methods of determining the relationship between cavity wall displacement and support pressure ("Ground response curve", GRC) for perfectly plastic and softening materials.Closed-form solution exist only for problems exhibiting rotational symmetry, that is circular tunnel-cross sections under plane strain conditions and spherical cavities.
Regarding the constitutive model, rigorous closed-form solutions exist for perfect plasticity, but only for a special case of softening, the Tresca material with linear strength decay 18 .This paper aims to fill this gap by providing a series of easy-to-use closed-form solutions for geomaterials that undergo strength reduction during shearing.The cases analysed are highlighted in bold in Table 1.Both cylindrical tunnels (represented by the circular hole model at plane strain) and spherical cavities are considered.Plane strain models, while simplifying reality by not accounting for the accurate stress-path followed by the ground during the advance of the tunnel face (see 20 ), are widely employed in tunnelling practice.This is because they allow for an intuitive and straightforward analysis of the ground-support interaction, known as the convergence-confinement method 21,22 , and Table 1.b) Neglecting completely the elastic strains in the plastic zone. (c) Salençon 15 . (d) Panet 16 . (e) Panet 16 , for the special case of zero residual cohesion. (f) Kovári 17 . (g) Carranza-Torres 18 . (i) This paper.

Perfectly plastic material
Softening material with linear decay law provide insights into the behaviour of tunnel and surrounding rock far behind the face.Spherical cavities are of importance as an approximation of the conditions in front of the tunnel heading 23 , providing insights into the development of ground deformation around the advancing tunnel face.In addition to the Tresca flow criterion (which is valuable as a model for clays under undrained conditions), the Mohr-Coulomb criterion is considered, which is more suitable for geomaterials under dry or drained conditions.The analysis of two strength decay laws (linear or exponential) allows a better representation of the observed softening behaviour.
The mathematical treatment can be significantly simplified by neglecting the elastic strains that develop during plastic yielding, i.e. assuming that the total strain inside the plastic zone is equal to the sum of the elastic strain at the elasto-plastic boundary and the plastic strain (Fig. 1a).Adopting this simplifying assumption (hereafter referred to as "Assumption A") reduces the mathematical complexity considerably (see 24 ), thus enabling closed form solutions.Panet 16 employed this simplification to derive closed-form Ground Response Curves (GRCs), for the specific case of a Mohr-Coulomb (MC) strain-softening material with zero residual strength.However, Guan et al. 25 conducted a comparative study on MC materials with constant friction angle and linear decay of cohesion and dilation angle in the strain softening regime, demonstrating that this simplification underestimates the radial cavity wall displacement by 20-40% compared to rigorous numerical methods (see 14 ).Furthermore, it can be proven that this simplification is incompatible with the zero volumetric strain constraint for materials following the Tresca yield criterion, in the case of spherical cavities, resulting in even more significant displacement underestimations (see Section "Error of the simplified solutions" of this paper).
Kovári 17 and Yu and Rowe 19 derived GRCs for perfectly-plastic MC materials by completely neglecting the elastic strains within the plastic zone (Fig. 1b).This simplification (hereafter referred to as "Assumption B") does not violate the incompressibility condition for Tresca materials, whether in relation to cylindrical or spherical cavities.However, it fails to ensure a continuous distribution of plastic strains and, assuming that strength in the strain-softening regime is controlled by the plastic strain, of strength, too.Specifically, disregarding completely the elastic strains in the plastic zone in combination with the condition of displacement continuity at the elastoplastic interface r = ρ results in a discontinuity in the plastic strain field; the plastic tangential strain at r = ρ + is by definition equal to zero while at r = ρ -it is equal to the elastic tangential strain at r = ρ + , that is ε t,ρ el .This discontinuity already results in a sudden drop in strength at r = ρ.
In this paper, starting with the definition of the problem (Section "Problem statement"), we adopt the aforementioned simplifying assumptions A and B (Fig. 1a and b, respectively) to derive in Section "Derivation" closed-form GRCs both for cylindrical and spherical cavities excavated in elasto-plastic strain-softening materials obeying the MC or Tresca yield criterion with exponential or linear strength decay (cases in bold in Table 1).Section "Synopsis of the ground response equations" summarizes and presents uniformly the derived equations.Subsequently, an interesting aspect of the predicted behaviour (the interplay between deformational behaviour and stability) and the effect of the softening law on the ground response are discussed (Section "Discussion").Finally, Section "Error of the simplified solutions" shows that the solutions derived adopting Kovári's simplification 17 (Fig. 1b) are, in most cases and particularly for spherical cavities, considerably more accurate than those derived under the simplification of Fig. 1a.

Problem statement
We consider a cylindrical or spherical cavity with radius a located in an infinite, homogeneous, and isotropic medium (Fig. 2).As in most GRC computations, the in-situ stress field is assumed to be isotropic and homogeneous.Solutions for anisotropic in-situ stress fields have been derived only for perfectly plastic materials [31][32][33] .
The assumption of isotropic in-situ stress is reasonable for deep tunnels 34 and sufficiently accurate for principal stress differences up to 30% 35 .
The cavity is supported by a uniform pressure σ a , while the far field stresses are equal to the in-situ stress σ 0 .Under these assumptions the system fulfils rotational symmetry.We are looking for the relationship between the support pressure σ a and the radial displacement u a of the cavity wall.The rock mass is taken as a linearly elastic, softening plastic material obeying the Hooke's law and the MC or Tresca yield criterion with linear or exponential strength decay (Fig. 3a and b, respectively) and constant friction and dilation angles.For the sake of simplicity, the possibility of out-of-plane plastic flow is not considered (see 27,31 and 32 for investigations into the effect of out-of-plane plastic flow in perfectly plastic materials and 1,2 for the case of softening materials).
Following 14 and 26 we assume that the decrease in strength is controlled by the deviatoric plastic strain γ pl , which is defined as follows: where ε pl t and ε pl r represent the tangential and radial plastic strains, respectively.In the case of linear strength decay (Fig. 3a): where f c , f p c and f r c denote the uniaxial compressive strength and its peak and residual values, respectively; δ is the maximum decrease in strength, normalized by the peak strength; and γ r is the deviatoric plastic strain at which the residual strength is attained.For exponential decay (Fig. 3b), Einav and Randolph's function 28 is adopted: where γ 95 is the deviatoric plastic strain at which the strength reduction equals 95% of the difference between peak and residual value.

Derivation
Upon cavity unloading, the rock mass remains initially elastic; the deviatoric stresses increase with decreasing support pressure and, at a certain support pressure σ ρ, may become so high that the stress path at the cavity wall reaches the yield surface and rock starts experiencing irreversible deformations.
The elasticity solutions for a spherical or cylindrical cavity in combination with the condition that the stress state (σ r , σ t ) at the cavity wall fulfils the yield condition with the peak strength parameters lead to the following expressions for the radial stress and for the strains at the tunnel wall at the onset of plastification (cf., e.g., 29 ): where m is a function of the friction angle (= 1 and (1 + sinφ)/(1-sinφ) for the Tresca and MC material, respectively); E and ν denote the Young's modulus and the Poisson's ratio, respectively; and ζ is the symmetry parameter (= 1 for cylindrical cavities and 2 for spherical cavities).
Further cavity unloading results to the formation of a zone of radius ρ where the rock mass experiences plastic shearing, and its strength decreases according to the softening laws of Fig. 3.In the case of linear strength (1) So, in the case of linear decay the system consists in general of three zones, hereafter referred to as the elastic zone (an infinite body with a cylindrical or spherical cavity), the softening and residual zone (thick-walled cylinders or spheres), and, in the case of exponential decay, the system consists of an outer elastic zone and an inner softening zone (Fig. 2).
The GRC can be determined by considering the aforementioned system components separately, taking account of the radial stress and radial displacement continuity conditions at their interfaces.
As mentioned in the Introduction, the simplifying Assumption A neglects the elastic strains that develop during plastic yielding, but takes account of the elastic strains that develop up to the onset of yielding (that is ), whereas the simplifying Assumption B considers the entire strain inside the plastic zone as being plastic, that is it completely disregards the elastic strains (ε t pl = ε t , ε r pl = ε r ).The derivation of the GRC will be illustrated considering the case of the MC yield criterion with linear strength decay and Assumption B. The solution for Assumption A can be derived analogously.

Kinematic considerations
Neglecting the elastic strains in the plastic zone allows the deformations to be analysed solely from the flow rule (ε r pl + ζ κ ε t pl = 0, where κ = (1 + sinψ)/(1 − sinψ) and ψ denotes the dilatancy angle) and the kinematic relationships, that is without making use of equilibrium and yield conditions 16,17,19 .
The assumption about the elastic strains in combination with the kinematic relationships and the plastic flow rule results in a differential equation for the radial displacement, the solution of which provides the displacement field in the plastic zone in the form u(r/ρ). So, in the case of Assumption B considered in this derivation, After we have determined the displacement field, we obtain from the kinematic relationships the radial and tangential strains and from Eq. ( 1) the plastic shear strain to be introduced in the softening law: Note that the form of the equations for the radial displacement of the wall and the distribution of the plastic deformation (Eqs.6 and 7) depend on which simplifying assumption is made (Eqs.6 and 7 hold for Assumption B), but not on the yield condition or the softening law.The latter are important for the relationship between radius ρ and support pressure, which is determined below based upon static considerations.

Static considerations
Introducing the yield condition where σ r and σ t denote the radial and the tangential stress, respectively, into the equilibrium condition with γ pl depending on r/ρ (Eq.7), the solution of which for the boundary condition σ r (a) = σ a provides the stress field in the form σ r = f(r/ρ, σ a ).This equation in combination with the condition σ r (ρ) = σ ρ results in a relationship between σ a and ρ, the specific form of which depends on the kinematic simplification (the aforementioned Assumption A or B), the yield criterion and the softening law.

Up to formation of the residual zone
Considering the initial (linear decay) part of the linear relation, the integration of Eq. ( 10) yields the following σ a (ρ)-relationship:

Onset of the residual zone
Equation ( 11) holds provided that there is no residual zone yet, which means that the plastic shear strain at the tunnel wall is less then γ r .This is the case (see Eq. 7) as long as the normalized radius of the plastic zone is less than a critical value ω: The support pressure at the onset of the formation of the residual zone (that is σ b ) can be obtained from Eq. ( 11) with ρ/a = ω.The introduction of the constraint "ω ≥ 1" in Eq. ( 12) is associated with the aforementioned discontinuity in the plastic strain field: If γ r is smaller than (1 + ζ κ) ε el t,ρ , then the plastic strain that is «instantaneously» imposed at r = ρ (Fig. 1b) suffices to cause formation of a residual zone right from the start, i.e. b = ρ.Τhat is for ρ/a > 1, and consequently ω (which according to the last r.h.s.term of Eq. 12 would become less than 1) must be enforced to the value 1.

After formation of the residual zone
If σ b > 0, then for σ a < σ b (or ρ/a > ω) a residual zone develops.The outer radius of the residual zone (b) is such that the support pressure exerted by the residual zone upon the softening zone equals to σ b (and consequently ρ/b = ω).It can be obtained from the bearing capacity equation of a fully plastified thickwalled cylinder or sphere that exhibits a uniform strength (the residual strength f c r ), is loaded at its outer boundary by σ b and is supported at its inner boundary by σ a (see, e.g., 17 ): The support pressure for ρ/a > ω is obtained from Eq. ( 13) and considering that ρ/b = ω, it reads as follows:

Closure
The GRC can be computed analytically by considering the displacement u a as an independent parameter, first computing ρ from Eq. ( 6) and finally σ a from Eqs. (11) or (14).
The derivation for the other cases (yield condition, softening law and assumptions about strains inside the plastic zone) follows the same path, the only difference being in the integration results.The next section summarizes the final equations in nondimensional form for all four combinations of yield function and softening law and the simplifying assumptions A and B. Vol

Synopsis of the ground response equations
The relationship between radius of the plastic zone and cavity wall displacement reads as follows: where ρ is the normalized radius of the plastic zone ( ρ = ρ/a ≥ 1 ); Ẽ and f p c are the normalized Young's modulus and peak uniaxial compressive strength, respectively, defined as follows: where the denominator is the in-situ stress (in the case of Tresca criterion) or its Caquot-transformation (in the case of MC yield criterion): The expressions for the support pressure depend on the specific combination of yield condition and softening law.They are given below using the normalization:

MC yield criterion, with exponential strength decay
where and E n () is the generalized exponential integral function.

Tresca yield criterion with exponential strength decay
where Ei() is the exponential integral function.

Discussion
Figure 4a and c illustrate the effect of the two softening parameters δ f p c and γ r Ẽ , presenting normalized GRCs for a bigger (upper diagram) and a smaller (lower diagram) strength decay ( δ f p c = 0.24 and 0.12 respectively) and a series of γ r Ẽ -values.The dimensionless parameter γ r Ẽ expresses how rapidly strength decreases with plastic deformations.The GRCs have been obtained considering MC yield criterion with linear strength decay, a cylindrical cavity (ζ = 1), φ = 25º, ψ = 5º and ν = 0.3, and making the simplifying Assumption B; comparisons with results for Assumption A and the rigorous solution are given in Section "Error of the simplified solutions".
The uppermost curve of every diagram applies to the limiting of a brittle plastic material characterized by a sudden decrease in strength as soon as the strain exceeds the strain at the peak strength.The lowermost line applies to the limiting case of a perfectly plastic material and is identical with Kovári's GRC 17 (a straight line in a double logarithmic diagram).With decreasing δ f p c -value (Fig. 4c), the effect of γ r Ẽ becomes predictably smaller as all GRCs move towards the GRC of the perfectly plastic material (lowest curve).
Figure 4b and d show the effects of the softening law on the GRC.The solid lines apply to the linear decay model, the dashed lines to the exponential model.To ensure a meaningful comparison between the two models, an equivalent γ 95 (exponential model) was calculated for every γ r (linear model).The estimate of γ 95 eq = 3/2 γ r results from the equalization of the integrals of the two functions, i.e. of the plastic works expended during shearing (shaded areas in Fig. 5).The exponential model with γ 95 eq shows a faster initial reduction in strength and therefore predicts larger deformations than the linear model.However, the difference is moderate and diminishes with decreasing δ f p c .As expected, the two models converge in the limit cases of perfect plasticity and brittle softening.
Due to the Caquot's stress-transformation, the GRCs extend to the range of negative support pressures which have no practical significance.Figure 4a and c show by means of vertical blue dashed lines the normalized support pressure σa at which the support pressure σ a becomes equal to zero for two cases concerning the post-peak strength decay: a partial strength loss (δ = 0.5) and a total strength loss (δ = 1).In both Figures, for δ = 0.5, all GRCs intersect the blue dashed lines, which means that the convergence of an unsupported opening stabilizes to a finite value, despite the decrease in strength.This is true for δ = 1, too, provided that γ r Ẽ is so big (e.g.48 and 120 in Fig. 4c) that a residual zone does not form.Otherwise, i.e. if a residual zone with zero uniaxial compressive strength develops, equilibrium is not possible 16 and, therefore, the displacements increase asymptotically to infinity as the support pressure tends to zero (e.g. the lines for γ r Ẽ ≤ 48 in Fig. 4c).The fact that in case of complete strength loss (δ = 1), equilibrium in the unsupported tunnel is not possible if a residual zone is formed, is evident from the GRCs of Fig. 6, where we refrain from Caquot's transformation in order to visualize better the behaviour at zero support pressure.The curves correspond to the normalized GRCs of the same colour of Fig. 4c, for specific values of σ 0 and E or γ r .The unsupported opening would be unstable in the case of the red solid line (δ = 1, γ r Ẽ = 24) and stable (i.e. in the case of smaller δ and/  or greater γ r Ẽ ).So, one may say that the normalized diagrams of Fig. 4 present in a unified way GRCs for cases that are qualitatively different from the constructional viewpoint (stable vs. unstable unsupported openings).
The stability condition of an unsupported opening in case potentially total strength loss (δ = 1) is σ b < 0, and can be determined in terms of γ r Ẽ and peak strength ( f p c ) based upon Eqs. ( 21) and ( 23), and it reads as follows: where Figure 7a shows this condition for varying values of φ (20°, 25°, 30°, and 35°) and ψ, where ψ is taken equal to φ-20°.It is interesting that an elasticity parameter, the Young's modulus, is decisive for the stability (or lack thereof) of an unsupported opening.This is because the Young's modulus affects the plastic strains (the latter increase linearly with 1/E in the case of a perfectly plastic material 30 ) and thus, in the case of softening, how much strength is lost.The opening type (cylindrical tunnel or spherical cavity) has a small effect on the stability condition.In Fig. 7b and c, the friction and dilation angles, respectively, are independently varied.An increased friction angle is advantageous for the stability of the unsupported opening, while on the contrary, an increased dilation results in larger plastic strains, leading to a more significant loss of strength in the case of a softening material.However, both φ and ψ have a relatively minor impact on stability, which justifies the simplification of keeping them constant in the strain softening regime.

Error of the simplified solutions
To assess the error of the simplified solutions, we employ the numerical model proposed by Lee and Pietruszczak 26 .This model enables the calculation of stress, strain, and displacement fields around a circular tunnel excavated in a strain-softening material.The plastic zone is discretized into a finite number of concentric rings and their thickness is determined to ensure equilibrium and compatibility.By using a sufficient number of annuli, this method accurately approximates the rigorous closed-form solutions.Figures 8 and 9 present the normalized convergence of an unsupported cavity for the cases of the MC and Tresca yield criteria, respectively.In both figures, the upper four diagrams hold for cylindrical cavities (ζ = 1) and the lower four diagrams for spherical cavities (ζ = 2).We investigate two values of the normalized peak strength, characterizing two squeezing intensities: moderate squeezing (f c p /σ 0 = 0.5, left diagrams) and heavy squeezing (f c p /σ 0 = 0.25, right diagrams).Two values for the maximum strength loss are considered: a smaller one, δ = 0.25, and a bigger one, δ = 0.75 for the MC material and 0.4 for the Tresca material.In the latter, a more unfavourable value of δ was not examined, since it would result in large strains and displacements that lack physical meaning and undermine the validity of the small-strain assumption, underlying both the analytical and the numerical methods.
Each diagram in Figs. 8 and 9 presents the dimensionless convergence as a function of Γ, i.e., γ r E σ 0 for the linear model and 2γ 95 eq E (3σ 0 ) for the exponential model.For each model, the following three curves are given; the rigorous solution based on the numerical method by Lee and Pietruszczak 26 (black curves) and the simplified solutions presented in this paper based on Assumption B (blue curves), and Assumption A (red curves), as derived from Eqs. ( 21), ( 22), ( 26) and ( 27). ( 32) . For the MC material (Fig. 8), both simplified solutions underestimate the convergence of the unsupported tunnel, with the proposed solution (blue curves) consistently improving Assumption A-based solutions (red curves).The improvement is particularly significant for the spherical cavity.Each curve approaches the solution for brittle plastic materials for small values of Γ and the solution for perfectly plastic materials for big values of Γ.The difference between the linear and the exponential model is insignificant and, as expected, visible only for intermediate values of Γ. Apparently, the numerical results (black curves) tend towards the exact borderline solutions, while the analytical results (red and blue curves) tend towards the simplified ones.
For the Tresca material (Fig. 9) and a cylindrical cavity (upper four diagrams), both simplified solutions closely approximate the exact one, with Assumption A slightly underestimating the convergence and Assumption B slightly overestimating it.Both simplified solutions converge towards the numerical one, as the analysis approaches the borderline perfectly-plastic and brittle-plastic cases, which are, in this case, unaffected by the simplifying assumption.Nevertheless, the difference is significant for the spherical cavity.While Assumption B only slightly overestimates the convergence, Assumption A greatly underestimates it.The latter is due to the violation of the zero volumetric condition, as discussed in the Introduction.
Figures 10 and 11 present the percentage error in radial displacement for the two simplified solutions, for MC and Tresca materials, respectively.
Despite the non-negligible errors present in both solutions, the proposed one (blue curves) consistently demonstrates reduced errors, remaining below 30% in most cases.An exemption occurs in the case of the cylindrical tunnel for the Tresca material, where the error of the proposed solution exceeds that of the Assumption A-based one.
However, unlike the Assumption A-based solution, the proposed one lies on the safe side.The error peaks for intermediate values of Γ, corresponding to the steep part of the curves depicted in Figs. 8 and 9.The maximum error is, in all cases, slightly smaller for the exponential model.eq E (3σ 0 ) for the exponential model, after 26 and according to the two simplified solutions (Tresca material). of qualitatively different cases, i.e. stable and unstable unsupported openings, in the same GRC.Specifically, for given cavity support pressure and strength drop δ, stable conditions are characterized by a finite dimensionless convergence of the cavity.Otherwise, the convergence tends towards infinity when the support pressure approaches zero.In the case of complete loss of strength, an unsupported tunnel cannot be stable.The stability condition for the latter case was presented in terms of peak uniaxial strength and Young's modulus normalized to σ 0 , and deviatoric plastic strain to reach the residual state.by Lee and Pietruszczak 26 .For the Tresca material and cylindrical cavities, both analytical solutions converge to the numerical one, with the proposed solution exhibiting slightly bigger errors, but on the safe side.A notable difference observed for spherical cavities, where the solution based on Assumption B slightly overestimates the response, while Assumption A significantly underestimates it.The linear and the exponential model, for equivalent values of γ r and γ 95 , predict similar convergences, but the exponential model exhibits a slightly smaller maximum error than the linear.

Figure 1 .
Figure 1.Radial distribution of tangential strain, elastic tangential strain and plastic tangential strain, (a), neglecting the elastic strains that develop during plastic yielding and, (b), neglecting completely the elastic strains inside the plastic zone.

Figure 2 .
Figure 2. Problem setup and zones of the stress field in a strain softening rock mass.

Figure 5 .
Figure 5. Linear vs exponential softening law for equivalent values of γ r and γ 95 .

Figure 9 .
Figure 9. Normalized convergence of an unsupported cavity over the dimensionless factor γ r E σ 0 for the linear model and 2γ 95eq E (3σ 0 ) for the exponential model, after 26 and according to the two simplified solutions (Tresca material).

Figure 10 .
Figure 10.Percentage error in radial displacement for the two simplified solutions (MC material).
Tangential elastic strain at the elasto-plastic boundary εel